
# Extract confidence intervals from parallel trend regressions
data_for_lwp <- graphs

# Set time around shock
data_for_lwp$time <- seq(-4,4,1)

# Creates points on each confidence intervals to search for LWP
Points <- matrix(ncol = nb_point , nrow = nrow(data_for_lwp))
for(i in 1:nrow(data_for_lwp)){
  Points[i,] <- quantile(c(data_for_lwp$lb[i],data_for_lwp$up[i]),seq(0,1,1/(nb_point-1)))
}

#############
# Search for coefficient from linear regressions using points on the intervals and verify whether fitted values are contained in the intervals

# Set parallel process
n_cores <- min(nb_point,detectCores() - 1)
n_cores
my.cluster <- parallel::makeCluster(
  n_cores, 
  type = "PSOCK"
)
registerDoParallel(cl = my.cluster)


# Regression 
round <- 0
Outs_all <- list()

Outs_all <- foreach(a = 1:nb_point, .combine = "rbind") %dopar%{
  
  Outs <- matrix(nrow = nb_point^8,ncol = 10)
  round <- 0
  
  for(b in 1:nb_point){
    for(c in 1:nb_point){
      for(d in 1:nb_point){
        for(e in 1:1){ #central point
          for(f in 1:nb_point){
            for(g in 1:nb_point){
              for(h in 1:nb_point){
                for(i in 1:nb_point){
                  round <- round +1
                  
                  nodes <- c(a,b,c,d,e,f,g,h,i)
                  
                  data_regression <- data.frame(time = seq(-4,4,1),points = diag(Points[,nodes]))
                  data_for_lwp$fitted <- lm(data = data_regression,formula = points ~ 0 + poly(time,dimension,raw = TRUE))$fitted.value
                  data_for_lwp$out <- ifelse(data_for_lwp$fitted <= data_for_lwp$up & data_for_lwp$fitted >= data_for_lwp$lb,0,1)
                  
                  Outs[round,] <- c(sum(data_for_lwp$out),nodes <- c(a,b,c,d,e,f,g,h,i))
                  
                }                  

              }                  

            }                  

          }                  

        }                  

      }                   

    }                  

  }                  
  Outs
}

stopCluster(cl = my.cluster)

# Keep outputs that fall within the confidence intervals
Outs_all <- Outs_all[!is.na(Outs_all[,1]),]
nodes_okay <- Outs_all[Outs_all[,1] ==0,c(2:10)]

# Keep the path with lowest coefficients on higher order term
magnitude <- vector(length = nrow(nodes_okay))
for(i in 1:nrow(nodes_okay)){
  data_regression <- data.frame(time = seq(-4,4,1),points = diag(Points[,nodes_okay[i,]]))
  R <- lm(data = data_regression,formula = points ~ 0 + poly(time,dimension,raw = TRUE))
  magnitude[i] <- R$coefficients[length(R$coefficients)]
}

# Run regression of LWP
data_regression <- data.frame(time = seq(-4,4,1),points = diag(Points[,nodes_okay[which(abs(magnitude) == min(abs(magnitude))),]]))

# longer and smoother path for graphing
R <- lm(data = data_regression,formula = points ~ 0 + poly(time,dimension,raw = TRUE))
lwp2 <- data.frame(path = predict(R,data.frame(time = seq(-5,5,0.25))),date = seq.Date(from = as.Date("2013-01-01"),to = as.Date("2023-01-01"),length.out = length(seq(-5,5,0.25))))

lwp2 <- lwp2[lwp2$date<= "2023-06-30",]
lwp2 <- lwp2[lwp2$date>= "2013-07-01",]


p <- ggplot(data=graphs, aes(x=date, y=coef)) + 
  geom_line(data = lwp2,aes(x=date, y=path),linewidth = 0.7,color = "#91959C",linetype=2)+
  geom_hline(yintercept=0, linetype=1, color = "black") + 
  geom_point(colour = "#B31B1B",size =3) + 
  geom_errorbar(aes(ymin=lb, ymax=up), colour = "#B31B1B", width=100,
                position=position_dodge(0.05))+ 
  labs(x="", y = expression(atop("Coefficients on the interaction between", "year FE and Change.Ded")),
       caption = paste("Polynomial order of the least wiggly path = ",dimension,"     -     ","pre-trend p−value = ",round(pt_value,3),"     -     post-trend p−value = ",round(pt_value_post,3),sep = " ")) + 
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  theme(legend.position = "none") 
